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Abstract 

The main purpose of the present paper is to study from a numerical analysis point 
of view some robust methods designed to cope with stiff (highly anisotropic) elliptic 
problems. The so-called asymptotic-preserving schemes studied in this paper are very 
efficient in dealing with a wide range of e-values, where 0 < e -C 1 is the stiffness 
parameter, responsible for the high anisotropy of the problem. In particular, these 
schemes are even able to capture the macroscopic properties of the system, as e tends 
towards zero, while the discretization parameters remain fixed. The objective of this 
work shall be to prove some e-independent convergence results for these numerical 
schemes and put hence some more rigor in the construction of such AP-methods. 

Keywords: Anisotropic elliptic problem, Asymptotic-Preserving scheme, Nu¬ 
merical analysis, Saddle-point problem, Inf-sup condition, Stabilization, Conver¬ 
gence. 


1 Introduction 

In a series of previous works [4, 5, 6, 11, 12] some efficient numerical schemes 
were introduced in the aim to solve at a moderate computational cost some highly 
anisotropic elliptic and parabolic problems. The interest in solving such problems 
comes for example from their regular occurrence in the modeling of magnetically 
confined plasmas [3, 9] and ionospheric plasmas [13], where the strong magnetic field 
creates anisotropy. An accurate and not resource demanding description of tokamak 
plasma dynamics is crucial for succeeding in the construction of a thermonuclear 
fusion reactor, producing clean energy for the future. 

The problems cited above involve a small parameter 0<£<1 measuring the 
anisotropy ratio in the diffusion matrix. This feature makes their numerical treatment 
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rather involved, since the problems degenerate in the limit e —>• 0 leading to a break¬ 
down of traditional schemes for e very small. This is caused both by the huge, 
e-dependent condition number of the discretized problem and by locking phenomena 
(the strong diffusion along a magnetic field line makes the solution to be almost 
constant along these lines, which is incompatible with an approximation by piecewise 
polynomials unless the computational mesh is well aligned with the field). In the 
previous works some efficient so-called asymptotic-preserving schemes were proposed 
and were shown to be able to cope with the deficiencies of traditional schemes. Their 
basic idea is to mimic on the discrete level the asymptotic behaviour of the continuous 
solution u £ in the limit e —>• 0, thus making the diagram in Fig. 1 commutative. 
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Figure 1: Properties of AP-schemes 


In the present paper, we restrict out attention to the case of elliptic linear prob¬ 
lems and are interested in two Asymptotic-Preserving schemes proposed in [6] and 
[12]. The efficiency and advantages of the different schemes was put into evidence 
numerically. However, the rigorous numerical analysis of these schemes is still lacking 
and is the subject of the present paper. The trick that makes these schemes work 
is the introduction of an auxiliary variable q £ which serves as a Lagrange multiplier 
in the limit e —>• 0 corresponding to the constraint on u £ , which results from the 
degeneracy of the governing equations. We are thus in the realm of mixed problems, 
their penalized variants and the discretizations thereof, as in [1, 8]. One cannot 
adapt though directly the techniques from these books to the present case as the 
inf-sup conditions are not satisfied on the discrete level, when one discretizes with 
standard finite elements as in the above cited papers. In fact, the choice of appro¬ 
priate functional spaces even on the continuous level is not straightforward. We are 
going here to propose an adequate functional setting with the inf-sup condition being 
introduced in a non standard way. We shall develop then a complete analysis of the 
finite element schemes with e-independent constants in the error estimates, relying 
on some discrete inf-sup conditions in /i-dependent norms. 

This paper is organized as follows. In Section 2, we introduce the anisotropic ellip¬ 
tic problem, which is the starting point of this work. A first Asymptotic-Preserving 
scheme for this problem, slightly modifying that proposed in [6] and well adapted for 
open field-line configurations, is then presented and analyzed in detail. Section 3 is 
concerned with the introduction of a different Asymptotic-Preserving reformulation 
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of the same anisotropic elliptic problem, being able to cope even with closed field¬ 
line configurations, which are often encountered in tokamak plasma modeling. This 
leads to the scheme proposed in [12]. A detailed numerical analysis of this scheme is 
then carried on. In Section 4 we validate numerically the error estimates obtained. 
Finally, some technical lemmas are postponed to the Appendices A and B. 


2 An AP-scheme for open field-line configura¬ 
tions 

Before presenting our model problem, let us first define some important quantities. 
Let b be a smooth field in a domain 17 C R rf , with d = 2,3, and let us decompose the 
regular boundary T = <917 into three components following the sign of the intersection 
with b: 

I'd ■= {a; € T / b(x) ■ n(x) = 0} , I N := I in U I out = {x £l / b(x) ■ n(x) ^ 0} . 

The vector n is here the unit outward normal to T. 

The direction of the anisotropy of our problem is defined by this vector field 
b € (C°°(H)) d , which is supposed to satisfy |6(a;)| = 1 for all x £ II. Given this 
vector field b, one can decompose now vectors v € M rf , gradients V0, with (f>(x) a 
scalar function, and divergences V • v, with v(x) a vector field, into a part parallel to 
the anisotropy direction and a part perpendicular to it. These parts are defined as 
follows: 

U|| := (v • b) b , v± := (Id — b® b)v , such that v = vm + v ±_, 

V|| (f> := (b ■ \7<£>) b, Vj_</> := (Id — b <g) b) V(p , such that V0 = + Vj_</>, 

V|| • v := V ■ ■un , Vj_ • v := V • v ±_, such that V ■ v = Vn • v + Vj_ • v . 

(1) 

Given these notations we can now introduce the highly anisotropic elliptic problem 
we are interested in, namely 

-iVy • (Ay V||tt e ) - V± • ( A ± V±u £ ) = f in Q, 

< \n\\ ■ (A||V||tt £ ) + nj_ • (Aj_Vj_rt £ ) = 0 onV N , (2) 

u £ = 0 on To. 

The parameter 0 < e < 1 is very small, inducing rather sever numerical difficul¬ 
ties, when solving (2) via standard methods. Indeed, this elliptic system becomes 
degenerate in the limit e —>• 0, leading to the reduced problem 

—V|| • (A||V||«) =0 in 17, 

n|| • (A||V|p) = 0 


(R) 


u = 0 
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on Tjv 
on r D 


(3) 




which has an infinite amount of solutions, all of them being constant along the field 
lines. Numerically this degeneracy translates in a very ill-conditioned linear system 
to be solved when 0 < e <C 1. 

The aim of the present section will be the mathematical study of the elliptic problem 
(2), in particular the investigation of its asymptotic behaviour as e tends towards zero, 
the introduction of an Asymptotic-Preserving reformulation, better suited to pass to 
the limit e —>• 0, and the detailed numerical analysis of the designed AP-scheme. 
The reformulation of the singularly-perturbed problem (2) is based on asymptotic 
arguments and is a sort of “reorganization” of the problem into a form, which al¬ 
lows for an automatic numerical transition from (2) towards the limit-model (to be 
determined) as e —>• 0, while keeping the discretization parameters fixed. 

2.1 Inflow Asymptotic-Preserving reformulation 

In order to avoid all the above mentioned difficulties corresponding to the non¬ 
uniqueness of the reduced problem (R), one has to pick up within all its solutions the 
right limit solution, by fixing in an adequate manner its value on the field lines. This 
was done in the previous works [4, 5, 6], via the introduction of Lagrange multipliers, 
which are necessary to recover the uniqueness in the limit e —>• 0. The numerical res¬ 
olution of the thus obtained Asymptotic-Preserving reformulations was shown to be 
stable and accurate independently on the parameter e, which is a great advantage as 
compared to standard discretizations for (2). This essential property of the designed 
AP-scheme was proved numerically, its rigorous numerical analysis being the subject 
of the present paper. 

For the mathematical study, let us assume that the diffusion coefficients and the 
source term satisfy the following hypothesis: 

Hypothesis A Let f £ iL~ 1 (P), 0 < e < 1 be a fixed arbitrary parameter and 


O 


T D ± 0. The diffusion coefficients Ay € VF 2,0O (P) and A± € M^ X rf(IT 2 ’ 00 (P)) are 
supposed to verify the bounds 


0 < Aq < Ay (x) < Ai, for a. a. x € Ll 


(4) 


Ao|M| 2 < v t A±(x)v < Ai||u|| 2 , Vu £ M. d with v ■ b{x) = 0 for a.a. x £ P, (5) 
with some constants 0 < Aq < A\. 

Before we shall pass to a brief presentation of an AP-reformulation of (2), we 
shall rewrite this problem in a slightly different form, masking the perpendicular 
derivatives, which turn out to be cumbersome for the numerical analysis. Indeed, 
the following reformulation, called in the following (P > ) £ -problem 




4jpV|| • (A||V||U £ ) - V- (AVu £ ) = / inH, 



(6) 


u £ = 0 


on T d . 


4 



is easily seen to be equivalent to problem (2) by setting 

A := (b <g) b) Ay (b <g) b) + (Id — b ® b) A± (Id — b®b). 

Remark that by Hypothesis A, we have immediately Ao||f|| 2 < v l A(x)v < Ai||r>|| 2 
for all v € and for a.a. x € which is a sort of coercivity and boundedness 
property for the diffusivity matrix A. 


Let us now introduce the mathematical framework and define the Hilbert space 
V as follows 

V := {v € , such that v\y d = 0}, (7) 

equipped with the scalar product 

(u, v)v = a(u, v) := / A'Vu-'Vvdx. (8) 

Jn 

In the following, the bracket (•, •) will stand for the standard L 2 -scalar product. We 
shall also frequently use the bilinear form an : V X V —>• M and the corresponding 
semi-norm | • L : V —> R defined by 



Ay Vim • Vuvdx , 


u 



(9) 


The weak formulation of problem (6) can be now written as: Find u e € V such that 

(P) e £ a\\(u £ ,v) + a(u e ,v ) = (f,v ), VT € V . (10) 

Thanks to Hypothesis A and to Lax-Milgram theorem, problem (10) admits a unique 
solution u e € V for all e > 0. 


The design of efficient schemes, which are uniformly stable along the transition 
e — >• 0, is based on the fundamental fact, that the solutions u £ € V of (10) tend for 
e —>• 0 towards some function n°, constant along the field lines of 6, i.e. belonging to 
the following Hilbert-space 

Q = {v € V , such that Vyn = 0} , (u, v)g := (Vj_u, Vj_u) L 2 (q) , (11) 

which consists of functions belonging to V with zero gradient along the field lines. 
Taking the test functions in (10) from Q , and passing formally to the limit £ —>• 0, 
permits to identify the problem satisfied by u° € Q, the so-called Limit model 

(L) a(u°,v) = (f,v) , Mv&Q. (12) 

Again, the Lax-Milgram theorem permits to show the existence and uniqueness of a 
solution u° € Q of this Limit problem (12). Remark that this Limit model is defined 
on a constrained space Q, and shall be equivalently reformulated in the sequel, on a 
constraint-less space. 
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The main idea behind the first AP-reformulation of problem (6) is to rescale the 
parallel derivative of u £ by introducing the auxiliary variable q £ such that 

V || q £ = ^V||U £ . (13) 

To ensure the uniqueness of q £ , we require in this section that q £ = 0 on Tj n . Remark 
that this is only possible if all the field lines are open and enter the domain (by Tj n ). 
For closed field lines, completely contained in 17, fixing q £ on Y, m would be not enough 
for the uniqueness and other methods shall be developed in Section 3. 

We thus introduce the Hilbert space 

Cin = {ge L\{ 7) / V|| q € L 2 (i 7) and q\ Tiri = 0} , (14) 

equipped with the scalar product a,|| (■, ■), inducing the norm | • y. Note that this is 
indeed a norm on Ci n since |g||| = 0 for q £ means Vy(/ = 0 on 17, which in 
combination with the boundary condition on implies q = 0. 

Substituting the definition (13) of q £ into (10) yields the following problem, called 
in the following Inflow Asymptotic-Preserving reformulation of ( P ) £ : Find ( u £ ,q £ ) € 
V x Cj. n satisfying 

f a(u £ ,v) + (1 - e)a\\(q £ ,v) = (f,v), VrG V 
(AP in ) £ \ (15) 

( tty (u £ , w) — ea» (q £ , w) = 0, Vtc € Ci n . 

The notation (APi n ) £ emphasizes the fact that we are introducing an Asymptotic- 
Preserving reformulation of (10), based on the Lagrange multiplier q £ , which is 
uniquely determined through the inflow boundary condition on Pj n . The reformula¬ 
tion (15) is completely equivalent to the starting model (10). However, remark that 
putting formally e = 0 in (15) and introducing for some technical reasons explained 
in the next subsection a larger space Ci n D Ci n leads to the well-posed problem : 
Find (u°,q°) € V x Ci n such that 

( a(u°,v) + a||(g 0 ,u) = (f,v), Vv € V 

(Lin) \ _ (16) 

[ ay (u^,w) = 0, Vw € Cin , 

which is an equivalent (saddle-point) reformulation of the Limit-problem (12). In¬ 
deed, instead of setting the problem on the constrained space Q, one introduces a 
Lagrange multiplier q° € C ln , enabling us to solve the problem on the constraint-free 
space V x 

2.2 The Inf-Sup condition 

Let us now focus on the mathematical study of the continuous problem (APi n ) £ and 
its asymptotic behaviour as e tends towards zero. Due to the saddle-point structure 
of this problem, we shall make use of the traditional inf-sup theory [1, 7]. The goal 
is to prove an e-independent inf-sup condition corresponding to (15), which ensures 
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the existence and uniqueness of a solution, as well as the convergence of the AP- 
solution ( u e ,q £ ) towards the L-solution (u°,g°) as e —> 0. For this, we shall need a 
more adequate norm on the space Anj in contrast to the one proposed in [6] (see (14)). 


Indeed, we would like that the form ay(-,-) satisfies an inf-sup estimate on the 
pair of spaces V plus the space of functions q. This would be trivially the case, if 
we would search for q in the space An> defined as the closure of An in the following 
norm |g|* 

a\\ (q,v) 

* := sup —j—:- 

uev Mv 

Note that for all q G C m , one has |g|* < |g| n , which means that the injection An C An 
is continuous, however An / An in general, as can be seen from the subsequent 
remarks. 



Remark 1 One can extend the continuous bilinear form to ay : An x V —>• M by 
defining for each q G An\An 

an(g, v) := lim a\\(q n ,v) , Vv G V , 

11 n— >■ oo 11 

for some {g n }neN C An such that q n ->n->oo q in An- 

Indeed, the sequence {g n }neN being a Cauchy-sequence in Ci n , one deduces immedi¬ 
ately that {ay (g n , n)} ne N Is also a Cauchy-sequence for each fixed v G V, being hence 
convergent. 

Remark 2 For any fixed q G F in , the maximum of a| j^j^ over v G V is attained 
with the function v* G V, which is solution to the problem 


(v* ,w)v = a\\(q,w) VwGV. 


(18) 


Indeed, let us fix q G Ci n and v* G V be the corresponding solution to (18). Then, 
any v G V can be decomposed as v = av* + v' with a := and v' G V verifying 

(n*,n')v = 0. We observe then that 


a\\{q,v) 

Mv 


a\v 


* 12 
IV 



V* 


12 
lv 


„/|2 


< 


\v |v = 


a\\(q,v*) 

|u*|v 


Remark 3 In order to prove that An An it suffices to verify that the norms in Ci n 
and Ci n are not equivalent, i.e. that there is no constant c > 0, such that |g|y < c|g|* 
for all q G C{ n . For this, it suffices to construct a sequence {qk}keN C An , such that 


a\\ {qk,v) ^ c . . c. 
sup — ,,, <7 =>■ \qk\* < T\qk 

vgV \qk\\\\v\v h h 


k 

Qk II ^ \Qk * j 
" c 


Vfc € N, 


(19) 


with c > 0 a constant. One can easily do it in the following simple setting: let 
hi = (0,7r) x (0, tt), An = 1, A |_ = Id, b = e 2 - Taking q^ = sin kx(cos y — cos 2 y), 
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which is a function in Ci n for any integer k > 0, we can find the solution v* k to 
problem (18) corresponding to q = qk as 

1 4 

v *k = sm kx cos y - cos 2 y). 

Now, in view of Remark 2 

a\\(qk,v) _ Q|| (qk,vj) _ \v* k \y _ J_ / 1 16 

veV l<?fc|||Mv |?fc||||u*|v \Qk\w Vb\k 2 +1 k 2 + 4' 

This gives an example of (19). 

Searching now for a solution ( u , q) belonging to V x £ in is the proper setting for 
our problem in the limit case e = 0. Indeed, in this particular case, we have to cope 
with a standard saddle point problem (Li n ) and the inf-sup condition is satisfied in 
the space V x Ci n - However, this choice does not work any more for e > 0, as the 
term ay(g, w) makes no more sense if we suppose only (q,w) € Ci n x Ci n . Hence, we 
propose to work for e > 0 in the previous space Ci n for the Lagrange multiplier q, 
associated however with the following slightly different norm 

Me : = (M* +e|g|jj)* , ( 20 ) 

which is equivalent to the old norm | • |y of Ct n with e-dependent equivalence constants 
exploding as e —>• 0 : 

Me < Vl + £|<7||| and My < 4= Me ■ 

v e 

The space (£j n , | • | e ) is a Hilbert one equipped with the scalar product 
(gi,92)e := iyl,v* 2 )v + £a\\(qi,q 2 ), Vgi, q 2 € C in , 

where v(, vl) are the unique solutions of problem (18). In the limit e —>• 0, this space 
transforms into the Hilbert space (£; n , | • |*) with the scalar product (qi,q 2 )* '■= 
{v\,v* 2 ) v . 


We are finally able to introduce the right mathematical setting for a rigorous 
study of the AP-problem (15) and its convergence towards the Limit-problem (16). 
The Hilbert space adapted to our problem is 


:= 


V X Cin for £ > 0 

V x Ci n for e = 0, 


u,q\\x e ■= (Mv + M* +eM|) 1/2 > 


and the problem we are interested in, can now be simply written as: Find for each 
e € [0,1] the solution [u £ ,q £ ) € X e to 


a{u £ , v) + (1 - er)a|| (<f, v) = (/, v), 
ay ( u e , w) — ean ( q £ , w) = 0, 

8 


M(y,w) € X e . 


(■ AP in ) £ 


( 21 ) 











For the further developments, we shall also introduce the coupled bilinear form C £ : 
X £ x X £ —> R defined as 


C £ ((u,q), (v,w)) := a(u,v) + (1 — e)a\\(q, v ) + a\\(u,w) - ea\\(q,w) . (22) 

This bilinear form C £ is uniformly continuous in e € [0,1], i.e. 

C e ((u,q), (v,w)) < 2\\u,q\\x e \\v,w\\ Xe , V(u, q), {v, w) G X £ , (23) 

as, using Cauchy-Schwarz inequality, one has 

C £ ((u,q),(v,w)) < |u|vMv + |<?|*Mv + MvM* +e|g||||io||| 

< ( 2 Mv + M* + e|9|j|) 2 ( 2 Mv + M* + eMjj) 2 • 

The form C £ enjoys furthermore the inf-sup property 


inf sup 

( u,q)ex e ( V:W ) £ x £ 


Ce({u,q ), (v,w)) 
\\u,q\\x e \\v,w\\x e 


>( 3 , 


(24) 


with a constant /3 > 0 that does not depend on e. This is established in the following 
lemma which is recast to a slightly more general and abstract setting. Our particular 
result is recovered from this lemma setting the bilinear forms a(-, •), b(-, •) and c(-, •) 
from the lemma to, respectively, a(-,-), a||(-,-) and ay(•,•). Note that this result is 
very close to those from Section 4.3 of [8] but we do not require here an inf-sup 
condition for the form b in V x L. 


Lemma 4 (Inf-Sup condition) Let V and L be Hilbert spaces with their respective 
scalar products o(-,-) and c(-,-) inducing the norms || ■ ||y and || • || l- Let moreover 
L D L be another Hilbert space with the norm || • ||^ such that \\q\\i < ||g||i for 
all q € L. Let b : L x V —» R be a bilinear form satisfying the continuity relation 
\\b(q, u)|| < IkHvIkllf, f or allv£V,q£L and 

inf sup = a > 0. (25) 

q&Lv&V IkllzJMIv 

Set L £ := L for any e > 0, Lq := L and let X £ for any £ > 0 denote the Hilbert 
space V x L e equipped with the norm ||u, q\\x e '■= (ll^lly + lkll|+ £ lklli) 1//2 - Introduce 
for any e € [0,1] the bilinear form C £ : X £ x X £ —>• R 

C £ ((u, q), (v, w)) = a(u, v) + (1 - e)b(q, v ) + b(w, u) - ec(q, w). 


Then C £ is continuous, with continuity constant M = 2, and satisfies moreover the 
inf-sup condition 


inf sup 

(u,q)GX e ( v , w ) £ x s 


C e ({u,q), (v,w)) 
\\u,q\\x e \\v,w\\ Xe 


>P, 


(26) 


with a constant fi > 0 that depends only on a. 
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Proof. To prove (26), let us fix an arbitrary (u,q) € X £ and denote 


Z : = sup 

(y,w)ex e 


C e ((u,q), (v,w)) 


\\v,w\\x e 

We want to prove that Z > /3\\u, q\\x e ■ First, we have 

/, \ I, II ^ (1^ C e {(u,q),(v, 0)) a(u,v ) ,, 

1 — e)a o t < sup-— - < sup-n - n-h sup -t— n— < Z + m v 

Mv v&v IMIv vev M v 


Now, we take v = u, w = — q and observe that 


C £ ((u,q),{u,-q)) = 


> 


a(u , u) — eb(q, u ) + ec(q, q) 

(i-|)IMIy + |lkll!, 


implying altogether 


1 

2 



+ 



ct 2 (l — e) 2 
8 



< Ce((u,q),(u,-q)) +^( Z + \\u\\ v ) 2 

< Z\\u, — q\\x e + + -|H|y- ■ 


Thus, 


1„ 

4 ^ 


2 + -I 
V + ol 


Q\\i + 


a 2 ( 1 — e)" 
8 


IMlf < Z\\u,q\\ Xe + -^ 2 < — \\u,q\\\ e + 


1 + 7 , 


for any 7 > 0 by Young inequality. Besides, we have for any e € [0,1] 


£ 

2 



ct 2 (l — e) 2 
8 


>co(lklli + e|klli)» 


with a constant Co > 0 depending only on a. Indeed, for e € [0,1/2] we can observe 
that (1 — £) 2 > | and conclude. For £ € [1/2,1], we can neglect the term with ||g||| 
on the left-hand side and use ||(?||^ < \\q\\L- 
Thus, assuming without loss of generality that cq < \ we have 


co\\u,q\\x s < — \\u,q 


+ i ± 2 Z 2 


Taking finally a sufficiently big 7 gives immediately ||«, q\\x e < (l//3)Z with a con¬ 
stant (5 > 0, independent of e. ■ 

The following theorem is the main theorem on the continuous level, which shows 
that the AP-reformulation (15) of problem (10) is well-posed and better adapted to 
capture the macro-scale behaviour of u £ in the limit e —>• 0. This AP-model provides 
thus a link between the micro-scale (e ~ 1) and the macro-scale (e ~ 0) behaviour of 
the system. 
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Theorem 5 (Existence/Uniqueness/e-Convergence) Let hypothesis A be sat¬ 
isfied. The AP-problem (21) is well-posed for each e G [0,1], i.e. for any f G V 7 and 
any e G [0,1] there exists a unique solution ( u £ ,q £ ) G X £ , which satisfies 

\\u/q e \\x e < 4ll/llv', 

with (3 > 0 the constant given by the inf-sup condition (26). Moreover, we have the 
e-convergences 

\\u e — u°, q £ — q°\\x e —>• 0 , for e —» 0 . 

If we suppose more regular data, as f G L 2 (fil), then one has even q° G Ci n and the 
estimates 

\u £ — -u°|y < Cy/e , \q £ - q°\* < C^/e, (27) 

with C > 0 some e-independent constant. 

Proof. The existence and uniqueness of a solution ( u £ ,q £ ) G X £ for each e > 0, is a 
simple consequence of the Banach-Necas-Babuska (hereafter BNB) theorem [7]. To 
prove the convergence ( u £ ,q £ ) —>• (u°,q°) we assume first that / G L 2 (Ll). We have 
proved in [6] that q° G Ci n in this case. Subtracting now (16) from (15) yields 

C £ ((u £ - u°, q £ - q°), ( v, w)) = ea\\ (q°, v + w ), V(w, w) G V x C in . 

Thus, for any e > 0, by the inf-sup property, there exist (v, w) G X £ = V x Ci n such 
that 

fi'\\u £ -u°,q £ - q°\\x e 11u, 11 x E < ea/q°,v + w) < e|g °|||\v + tn||| < v / 2e|g 0 ||| ||n, w\\ Xe , 

with some 0 < fi' < ft, for ex. ft' := fij2 , which implies \\u £ —vP, q £ — q°\\x e < k°|||, 

leading to the convergence estimates (27). 

We are now going to generalize this result to any / G V' by a density argument. Let 
us denote simply by U £ (f ) the solution ( u £ ,q £ ) G X £ of (21) associated to / G V 7 . 
Now fix some / G V'. Since L 2 (fiX) is dense in V 7 , for any 5 > 0 there exists fs G L 2 (Ll) 
such that f = fs + Rs with ||i? 5 ||v' < Hence, there exists £o > 0 such that for 
all e < £o 

\\u e (f)-u°(f)\\ Xe < \\u £ (fs)-u o (fs)\\ XE +\\u £ m-u°m\\ Xe < 6 -+j\\r5\\v< 5. 

Here we used the fact that fs G L 2 {fil) which implies || U e (fs) — U°(fs)\\x e < —> 0 

as £ —>• 0. ■ 

2.3 The numerical analysis of the inflow AP-scheme 

Having reformulated on the continuous level the singularly-perturbed problem ( P) £ 
into a system (APi n ) £ which is better suited to capture the macroscopic limit as 
£ —^ 0, we shall now discretize via a standard approach this new system and analyse 
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the obtained AP-scheme in detail. In particular error estimates are deduced and the 
convergence of the scheme independently on the anisotropy parameter e is shown. 

Let us introduce a mesh 74 on hi consisting of triangles (resp. rectangles) of 
maximal size h, let 14 C V be the finite dimensional space of P& (resp. Q&) finite 
elements on 74, and let us define Lh := 14 H Ci n = 14 , n £,; n as well as Xh := 14 x L^. 
Note that we require Lh C 14, which signifies that we enforce the boundary conditions 
on Tp for functions in Lh, cf. [6]. We are thus looking for a discrete solution 
{u £ h ,q £ h ) € V h X L h of 


(4Pj n )4 


a{u £ h , v h ) + (1 - e)a\\ {q £ h , v h ) = (/, v h ), Mv h <E 14 
ay [u £ h , w h ) - eay (q £ h , w h ) = 0, \/w h € L h . 


(28) 


The analysis of this scheme would be straightforward if the discrete inf-sup condition 


. . a\\(qh,Vh) 

mt sup -———— 

q h &L hVhe y h \qh\*\Vh\V 


> a , 


(29) 


were satisfied with an e- as well as mesh-independent constant a > 0. However, 
this constant is unfortunately mesh dependent, as shown in Appendix C. In order to 
circumvent this difficulty we introduce the following mesh-dependent norm on L} t 

a\\(qh,v h ) 

*h • Slip : ; • 

v h ev h \Vh\v 


Note that this is indeed a norm, since \qh\*h = 0 implies an (<&, (fa) = 0 due to the 
inclusion Lh C 14 and thus Vn(fa = 0, which in combination with the boundary 
conditions on r* n yields qn = 0. 


We now equip the space Xh with the norm 

IK, qh\\x S ' h = (\u h \v + \qh\l h + e\qh \\) lf2 ■ 

By Lemma 4, the bilinear form C e is continuous on Xh with this norm, and enjoys 
the inf-sup property 


C e ((u h ,qh ), (vh,w h )) 

( u htQh)^^h (vh ,Wh)£Xh W'U'h, Qh II 11^, l|x £i j 


inf 


>p, 


(30) 


with a constant /? > 0 that does not depend neither on the mesh nor on e. This 
implies the discrete version of theorem 5. 


Theorem 6 (Discrete Existence/Uniqueness/e-Convergence) The discrete AP- 
problem (28) admits for each fixed h > 0 and e > 0 a unique solution {u £ h ,q £ h ) € 

V h x L h , satisfying 

KKlk.h < ^II/IIv', 

and one has the e-convergence 

II4 - K Qh - Qh\\x s , h ^0 for e -> 0. 
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Moreover, the condition number of the matrix corresponding to problem (28) is 
bounded by a constant independent of s (assuming that the same bases of Vh and Lh 
are chosen for all values of e). 


Proof. The existence and uniqueness of a solution ( u e h , q e h ) € Xh for each e > 0 , is a 
simple consequence of the BNB theorem [7]. The convergence kk) —kk) 
can be established by the same arguments as in the proof of Theorem 5. 

We turn now to the study of the condition number. Let {(f)f, .. ., (resp. 

{4 >\,..., f g Nq }) be a basis of Vh (resp. Lh). We shall identify every function Uh € Vh 
(resp. qh € Lh) with a vector u € H Nu (resp. q € M Nq ) consisting of the expansion 
coefficients of Uh (resp. qh) in these bases. Denoting the Euclidean norm of a vector 
by || • 112 and using the equivalence of norms on a finite dimensional space, we observe 
that for all Uh € Vh and qh € Lh we have 


hu\\u\\2 < \Uh\v < DilMk , hq\\q\\2 < M|| < Vq\\q || 2 , 11 ^| 1 2 < Wh\*h < V*M\2 


with some positive constants /r’s and z/’s. We shall moreover identify any = 
( Uh,qh) € X h with a vector € R N , N = N u + N q , such that = (u r ,q T ) T . We 
observe for any such < 1 ^ that 11 I 2 = IMIl “HMl!) which in combination with the 
estimates above gives 


min{/i u ,/z* +e ^}||$|| 2 < ||Mk> < niax{z/ u , 14 + ev q 


2l "$lli 


Let now A denote the N u x N u matrix with entries = a(0“, </>“), B the N u x N q 
matrix with entries bij = , (fj), and C the N q x N q matrix with entries Cij = 

a\\(cj) g , (j) g ). The matrix corresponding to problem (28) can be then written in the 
following block form 

A (1 -e)B 
B t eC 


A e = 


Its 2-norm denoted by 11 • 11 2 is bounded for all e € [0,1] by 

...... ’L-A £ <3> C £ (^ h ,^h) 


2 - 


sup 


ll'hlUI'hllz 


sup 


< M sup 




|<L>| | 2 I 


*h,*heXh\{ 0 } ||^|| 2 ll^ll 2 

e ' h < M rnax(^, v) + v q ) 


where M is the (e-independent) continuity constant of C e . Similarly, using the inf- 
sup property of this bilinear form we know that for all € Xh there exists T/, € Xh 
such that 

^mm{/4,^}||$|| 2 ||tf|| 2 < kMkJIMk,* < C e (Q h , *„) 

= T • A e <L < | k11 2 11A e ^11 2 , 

with an e-independent constant (3 > /3 1 > 0. This simplifies to 

/3'nl}\\$\\ 2 < ||A e $|| 2 , 
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or equivalently 


/3'i4}\\(A £ ) 1 $\\ 2 < ||$|| 2 , Vlel^ 
Thus, the condition number can be estimated as 


cond 2 ( A £ ) = ||A e || 2 ||(A l 


£\ — 1 I 


|2 < 


M max(i/‘, v( + v 2 ) 
P min{/i2 , p 2 } 


which is an e-independent bound. 


( 31 ) 


Remark 7 Let us try to be more quantitative in our estimate of cond 2 (A £ ). In 
what follows, the symbols < and ~ will hide the constants of order 1, independent of 
the mesh. Consider the standard finite element setting: the bases of \f and L^ are 
formed by the hat finite element functions on a quasi-uniform mesh. We know in 
this case that ||u/i ||| 2 ~ 11 ^ 11 i an d \v,h\v < Cjh^ 1 ||tt/i||z ,2 by the inverse inequality 

with a constant Cj > 0 that depends only on the mesh regularity [7]. We also recall 
the Poincare inequality ||rtfe || L 2 < Cp\uh\v- The same holds for qy, and leads to 


hu ~ hg ~ h an d v u ~ v q ~ 1 • 


We also have \qh.\*h. < M|i, hence v* < v q . Moreover, for any qh € Lh we prove, 
using the inverse and Poincare inequalities, that 




a\\(qh,qh ) 


Ml 


\qh\v 


(Mj_ + IMjj) ' (Cih- 2 \\q h \\ 2 L 2 + Mjj) 


> 


Mi 


1/2 


C 2 P C 2 jh 2 \q h \l + Mi 


> 


(C 2 C 2 h~ 2 + 1) ( C 2 h -2 
{C 2 h~ 2 + C 2) 1/2 


2 

L 2 


1/2 


> 


C]h 2 \\q h \\ 2 L 2 + \qh\ 


1/2 


C 2 P C 2 h~ 2 + 1 


Ml 

MIl 2 ~ h\\q h \\ L 2 ~ /i 2 ||g || 2 . 


C 2 Cjh~ 2 + 1 


This implies /i* > h 2 , so that (31) becomes finally 


cond 2 (A £ ) < 


Theorem 8 (h-Convergence) Let k > 1 and Vy, C V be the P& (or Q^) finite 
element space on a regular mesh Th- Suppose moreover that problem (21) has a 
solution ( u £ ,q £ ) € X £ , having the regularity u £ € H k+1 (Ll), q £ € H k+1 {Ll). Then one 
has the estimate 

\u £ - u £ h |v < ch k (\u e \ H k+i + |<fM+i) i (32) 

with a constant c > 0 that depends neither on the mesh, nor on e. 
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Proof. Let u £ h € V/,, and q £ h € be the standard nodal interpolant of u £ and q £ 
satisfying [7] 

| u £ — u £ h \ H i < ch k \u £ \ H k+i and \q £ — q^lm < ch k \q e \ H k+i . 


We can now derive the error estimates in the H 1 -norm for u £ in the way similar to 
Cea’s lemma: by the inf-sup property, there exists ( Vh , w>h) € Xh with \\vh, Wh\\x e h = 
1 such that (with some 0 < /3' < /3) 


u h ~ Uh\m 


< < jyC £ ((ul-u £ h ,q £ h -q £ h ),(v h ,w h ))(33) 

= ^C £ {{u £ -u £ h ,q £ - q £ h ),(v h ,w h )) 

< c(\u £ - u £ h \l + | q £ - q £ h \l h + e\q £ - q £ h |J) 1/2 

< c(\u £ — u £ h \ H i + \q £ — qllni), 


since 


, a\\(q £ -q £ h ,v h ) 

\*h = sup - i—i- < I q - " 

v h ev h \Vh\v 


< \q £ ~ 


I m- 


We can now employ the interpolation error estimates to conclude. 


Remark 9 The error estimate (32) would be of course useless if the norms \u £ \ H k+i, 
\q £ \ H k+i were e—dependent and exploding in the limit £ —> 0. Fortunately, it is not 
the case. We expect indeed that \u e \ H k+i is bounded uniformly in e by the norm of 
f in H k ~ 1 (Q) and \q e \ H k+i is bounded uniformly in e by the norm of f in H k+1 (Q). 
This can be easily proved in the case of a simple aligned geometry, see Appendix A. 
We conjecture that this remains true also in a general setting. 


Remark 10 If we do not omit the norm of q ^ — q^ in the left-hand side of the first 
inequality in (33), we also get an error estimate for q £ h 


h k 

\q £ — q.h\\\ < c-j=(\u e \ H k+i + \q £ 




which degenerates as e goes to 0. We are not sure, if this estimate is sharp, but we 
recall that q £ is an auxiliary variable, without any intrinsic meaning. 


3 Second AP-reformulation for general field lines 

The fundamental idea of the AP-reformulation introduced in Section 2 is the intro¬ 
duction of a Lagrange multiplier q £ E Ci n in order to handle well with the constraint 
V||it° = 0 in the limit e —» 0. This Lagrange multiplier was uniquely determined up 
to a constant on the field lines, which was fixed by imposing q( r = 0. The disad- 
vantage of this scheme is that it requires to identify the inflow part of the boundary, 
which can be cumbersome in practice or even not possible if some of the field lines 
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are closed and lie completely inside the domain Q. It is thus tempting to abandon 
the zero inflow boundary condition and to search for the auxiliary g £ -variable in the 
Hilbert space 

£ = ^eL 2 (H) / V||£€L 2 (H)}, (u,v) c :=(u,w) + (V\\u,V\\w). ( 34 ) 

The problem with this idea is that we loose now uniqueness of the solution if we 
attempt to implement the AP-reformulation (15) just changing to C. To cir¬ 
cumvent this difficulty, it was proposed in [ 12 ] to introduce a stabilization term into 
the AP reformulation so that it becomes: Find (u £,a , £ e,fT ) € V x C such that 


(AP S ) 


£,cr 


a(u £,<T , v) + (1 — e)a\\ (£ e ’ CT , v) = (/, v), Vu € V 

oil ( u £ ’ a , w) — ea\\ (£ e,cr , w) — cr{i £ ' a , w) = 0, \/w € C , 


(35) 


where cr > 0 is a small stabilization parameter, chosen consistently with the overall 
discretization error. It is this term which permits to have the uniqueness, as will be 
shown in Lemma 13. 

In the limit e —>• 0 this system yields: Given a > 0 , find (u 0,cr , £ 0,fT ) € V x £ 2 , 
solution to 


{LsY 


a{iY a ,v) + ay (Y ,rT ,v) = (/, v), Vv € V 
oil (u 0,rT , w) — fj(^ 0,cr , w) = 0, \/w € C 2 , 


(36) 


where C 2 is, loosely speaking, the closure of C in the | • |* semi-norm (17) intersected 
with L 2 (H), i.e. 

C 2 = < ^ € L 2 (H) / sup 
l i)GV 

This space is a Hilbert-space associated with the scalar product 


Q||(g^) 

Mv 


< oo 


(u, U ))£2 := (u, w) + (u*,w *) L 2 , V(u, w) € C 2 , 


where u* resp. w* are the unique solutions of (18) corresponding to u resp. w. We 
need this special space, first of all, to be able to treat the limit-problem (Ls) a with 
the inf-sup theory, similar to the inflow-case, and also in order to be able to define 
the stabilization term cr(£ 0,a ,w). 


Remark 11 Remark also that we have C 2 ^ C. Let us prove it in the following 
simple setting: let Ll = (0,7r) x (0,7t), Am = 1, A± = Id, b = ei- For any q = 
YYki =l Qki si n koc cos ly, the calculation as in Remark 3 gives 


W\ 


2 

* 


E 


k,l =1 


Z 4 

k 2 + Z 2 qkl 


2 


so that taking qki such that qki 


7 if k = l 2 and q k i 


0 for any k / l 2 we have 


00 72 

\ q \ 2 * = J2 H ^]2 

1=1 


< oo , 


16 





so that q G C. Moreover, clearly q G L 2 (Gl). However, 

OO 

/ 2 M 2 = °o. 

fc,Z=l 


3.1 Mathematical analysis on the continuous level 

To analyze the well-posedness of problem (35) ans its asymptotic limit behaviour for 
e —x 0 , we shall rewrite it in an equivalent manner, better suited for mathematical 
studies. For this, we observe first that the second equation in (35) gives for all e, a > 0 

(£ £ ’ a ,w) = 0 , VweGc, with Gc ■= {v £ C \ Vp = 0 } , 

which means that £ £,cr belongs to the following space 


A = {£ e £ / (£, w) = 0 \/w g Gc ), 


(37) 


which consists thus of functions from C with zero (weighted) average along the field 
lines. Remark that one has the decomposition C = Gc. ©"*“ A. Problem (35) can be 
hence rewritten as: Find » £ £, °") G V x A such that 


f a{u £ ' a ,v) + {l - e)a\\(£ £ ' a ,v) = (f,v), Vu G V 

{AP's) £ '° \ ( 38 ) 

[ a\\(u £ ' a ,w) — ea\\(^ £,rT ,w) — a(£ £ ’ a ,w) = 0, VwgA. 

We emphasize that this reformulation is completely equivalent to (35) for all e > 0 and 
a > 0 and is done solely for the purposes of mathematical analysis. The formulation 
used for the numerical discretization will be (35). 

Note that A becomes a Hilbert space when equipped with the scalar product 
<X|| (-, •) and corresponding norm | • y. Indeed, if £ G A and |£||| = 0 then £ G Gc which 
implies £ = 0 since A is orthogonal to Gc- For the same reasons, the semi-norm | • |* 
(17) is actually a norm when applied to A. We can thus introduce the closure A of A 
with respect to | • |*, needed as usual, for the e —>• 0 limit model. The (Ls^-problem 
will be shown to be equivalent to: Find (u°’ a , £ 0,<T ) G Vx (in L 2 ( fl)), solution to 

( a(u°’ a ,v) + a^°’ <T ,v) = Vu G V 

{L'sY { _ ( 39 ) 

| ay (u°’ a , w) — w) = 0, Vu; G A D L 2 (Gl). 

As mentioned earlier, formulations (38) and (39) are better adapted for the math¬ 
ematical study, then the completely equivalent ones (35) and (36). In particular in 
the limit er —x 0, they permit to get the following problems: Find (u e ,£ £ ) G V x A 
solution of 


f a{u £ , v) + (1 - £>11 (£ £ , v) = (/, v), W G V 
(AP a ) £ { (40) 

[ ay (u £ ,w) — £d||(£ £ , w) = 0, Vw£A. 

which is equivalent to the original problem (10) and hence also to the inflow AP- 
reformulation (15). In the present case, we fix the Lagrangian variable £ £ by requir¬ 
ing zero average along the held lines, i.e. G A, in the former case we fixed the 
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corresponding Lagrangian variable q £ by setting q £ zero on the inflow boundary L,; n . 
i.e. q £ E Ci n . Note that we do not want here to discretize the space A directly. This 
space arises only in the limit a —>• 0 , which is never performed in practice when one 
implements the scheme of this paper. On the contrary, the scheme from [ ] relies on 
a direct discretization of A which results in a rather complicated numerical method. 
Remark also that we abandoned in (40) the requirement that the ^-variable has to 
belong to L 2 (fl), as there is no more need, for a = 0 . 


Letting now formally e —> 0 in (40), we obtain the problem: Find (u°, £°) E V x A 
such that 

( a{u°,v) + a||(£°,u) = (f,v), \/v E V 

(La) { n - (41) 

( ay (u u , w) = 0, Vu> € A , 

which is an equivalent (saddle-point) reformulation of the original limit problem ( 12 ). 


For the reader convenience, we draw in Figure 2 a scheme, with all the problems 
we introduced so far, and their relations. In the following Lemmata and Theorems, 
we shall prove some of these relations and convergences, adapting the results from 
the previous section 2.2 to the present case containing two parameters, e and a. 


equiv. 


o —y 0 


equiv. 


(py 


0 ap a y 


AP' s ) e ' a 


C AP S Y'° 


£ —y 0 

£ —y 0 

£ —y 0 


o —y 0 

equiv. 


{La) 


{L's)° 


{Ls)° 


Figure 2 : Stabilized reformulations of the original problem (P) £ . 


Lemma 12 (Inf-Sup condition) Let V, L, L, L be Hilbert spaces such that L C L 
and L C L with continuous inclusions and ||^||^ < ||£||l for all £ E L. Let a(-,-), 
c(-, ■), d(-, •) denote the scalar products on respectively V, L, L and b(-, ■) :LxL ->R 
be a bilinear form satisfying ||6(^,u)|| < ||f||y||C||^ for all v E V, £ E L as well as the 
inf-sup condition 


mf sup —. 

ZeLvev IrIIlFIIv 


= a > 0 . 


( 42 ) 


Define furthermore the Hilbert space X Ey(J for e > 0, a > 0 by 


X E 


V x L, if e > 0, <7 > 0 

< V x (L n L), if e = 0, o > 0 

V x L, if e = 0, o = 0 
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and equip it with the norm ||it,£||x £>CT := (IMIy + ||f ||| + e\\(\\ 2 L + <j\\i\\)) 1 / 2 ■ 

For any e > 0 and a > 0 let C E , a : X £A x X St(T —> M be the bilinear form defined 
by 


C e ,a((u, £), (v, w)) = a(u , v) + (1 - e)b(£, v ) + b(w, u) - ec(£, w) - ad(£, w). 
Then C £j(7 is continuous and satisfies the inf-sup condition 


inf sup 

(u,OS-Y e , CT (v,w)eX s , a 


CeA( U iOi (v,w)) 

\ u ^\\x e A\v,w\\x^ 


>p, 


with a constant /3 > 0 that depends only on a. 


(43) 


Proof. The proof of this lemma follows the same lines as that of Lemma 4 and we 
give here only a short version of it. For any (u,£) E X Et(T , denoting 


Z : = sup 

(»,iu)6X e , 5 


\\v,w\\ Xei<7 


we can prove that (1 — e)a||^||^ < Z + ||it||y. Now, taking v = u, w = — £ we observe 
that 


CeA( u >0>( u >~$)) = a (u,u) - eb(£,u)+ £c{f,,f,) + crd(£,£) 

> ( 1 - £ AMlr + £ M 2 L + °ml, 


implying altogether 

i|MII + f Kill + ° 2(1 ~ 6)2 Kill + allflll < Z||»,-«ilx„ + A + . 

Following again the inequalities from the proof of Lemma 4, we see that there exists 
a constant cq E (0, |] depending only on a such that for any 7 > 0 and e E [0,1] 


colkClIx,,. < T||mKII 


I +i ± L? 2 . 

XV e u O 


Taking finally a sufficiently big 7 yields ||u,£||x s < (1 /fi)Z with a constant f3 > 0 
depending only on a. M 


Lemma 13 (Existence/Uniqueness for e > 0 and a > 0) Let hypothesis A be 
satisfied. The stabilized AP-problem (APs) £,a (resp. (Ls) a ) is well-posed for each 
e E (0,1] and a > 0 (resp. e = 0,a > Q), i.e. for any f E V there exists a unique 
solution (u £,a , A’ a ) E V x C (resp. (u 0 ’ a ,(°’ a ) € V x C 2 ), which satisfies 


\u 


£,<J 


,£^11* <C||/|| V , 


(44) 


with ||tL£lk, CT --=(Mv + l£l* + e l£ly + a U \\ 2 L 2 ) 1 / 2 and some C > 0 independent on e 
and a. Moreover we have the e- convergence 

\\u £ ’ a -u°’ a ,A’ a ~ A ,rr \\x e:!T ->0 for e -> 0 . 
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Proof. The existence and uniqueness of the solution to the reformulated problems 
(. AP' s ) £,a and ( L' s ) a follows directly from Lemma 12 by setting V = V, L = A, 

L = A, L = L 2 {fl). Now, the equivalence of ( APs) £,a and (AP' s ) £,a is easily seen 
from the decomposition C = Qc A. Similarly, the equivalence of ( Ls) a and ( L' S Y 
can be derived from the decomposition C = Qc ©^ A. ■ 

Theorem 14 (Existence/Uniqueness for e > 0 and a = 0) Let hypothesis A be 
satisfied. The (APj/) £ -problem (40) (resp. ( L/)-problem (41)) is well-posed for each 
e E (0.1] (resp. e = Q), i.e. for any f E V and any e E [0,1] there exists a unique 
solution ( u £ ,f £ ) E which satisfies 

ll u£ >£ £ IU e ,o < C'll/llv'> 

with some C > 0 independent on e. Furthermore, one has the e-convergence 
IK* £ - *i°,£ £ - f°||* e ,o -4 0, for £ —>• 0 . 

If we suppose more regidar data, as f E L 2 (Q), then one has even £° € A and the 
estimates 

\u £ - rt°|y < C/e , |£ £ - £°|* < C/e , 

with C > 0 some e-independent constant. 

Proof. The existence and uniqueness of a solution (rt £ ,£ £ ) E V x A to (APX) £ resp. 
(ii°,£°) € V x A to (Ljfi) is easily established using Lemma 4. The statements about 
the convergence as e —>• 0 follow in the same way as in the proof of Theorem 5. ■ 

Theorem 15 (a-Convergence) Let hypothesis A be satisfied and moreover, (u £,a , £ £,<T ) E 
V x A be solution to ( APs) £ ’ a and (it £ ,£ £ ) E V x A solution of [APj/) £ , with e > 0. 
Suppose that £ £ E Then 

\\u E ’ a -u e ,e’ a -e\\x s , a <ca\e\m, ( 45 ) 

with a constant c > 0 independent of a and e. 

Proof. We turn now to the convergence as cr —» 0. Using the combined bilinear form 
C £jtT and recalling the problems ( AP' s ) £,a resp. (AP_a) £ , we can write 

C £ia ((u e ’ a , £ £ ’ a ), (v, w)) = {f,v), V (v, w) € V x A , 

C £ ,o((u £ ,£ £ ),(v,w)) = (/,«), V(w,ra)eVxd, 

where 

C e , a {{u,£), (v,w)) := a(u,v) + (1 - e)a||(£,u) + a/u,w) - ea/(,w) - cr(/w). 

Note that C £) o coincides with C £ as defined by (22). Taking the difference gives 
C £ , a ((u e ’ a - u £ ,( £ ' a - £ £ ), (v,w)) = <r(£ £ , w) < a|£ £ |y|u;|y/ < c<r|£ £ |y|u;|*. 
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We have used here the bound |w|v' < c|w|* valid for w € A as proved below (Corollary 
19). 

Now, remind that the form C £iCT enjoys the inf-sup property 


inf 


sup g^( ( ". 0 . (l, .« ,)) 


X II /-II II II 

(w,^)GA’e )£T (v,w)(iX II^J S ||<V e> o- 11^’ ^ll^e 

where ||u,£||* ei0 . = [\u\$ + |£| 2 + e|£|g + cr||C||^ 2 ) 1/2 , so that H* < \\v,w\\ Xe , a - We 
can thus conclude that there exists (v,w) € X £)CT such that ||u, w\\x S:t7 = 1 and 

/3'\\u £ ’ a -u £ A £ ^-C\\x e , tT < C £ta ({u £ ’ a -u £ ,C £ ’ a -^ £ ),(v,w)) < ca\f, £ \v\\v,w\\ Xe rr = ccr|£ £ |^i 

with some 0 < f}' < (3, for example /3' = (3/ 2. This concludes the proof. ■ 


Remark 16 Without the additional hypothesis f £ G i7 1 (0), we can easily prove a 
sub-optimal estimate 


\u 


e,a 


\\x. a < CyftU e \\ L 2 


Indeed, 


C £A ({u £ ' a -u £ ,£ £ ' a -i £ ),(v,w)) = a(f £ ,w) < CT||^ £ || I; 2 (n) ||u;|| L 2 (Q) < cy/a\\^\\ L 2^\\v,w\\ Xe , a ■ 

Remark 17 The conclusions of Theorem If remain true (after an obvious rephras¬ 
ing) in the limit case e = 0 since the proof relies on the estimates in the norm of X £)(T 
which remains a valid norm in the limit e —>• 0. 


It remains to prove the bound |rr|y' < c|u;|* valid for w G A. This will be done 
using the following result: 

Lemma 18 Let u G V and consider v & A being the unique solution to 

a\\(v,w) = (u,w), Vic G A. (46) 

Then v G H 1 (fl) and there exists a constant c > 0 such that ||u||/p < c|tt|v- 

Proof. To simplify the notations, let us restrict ourselves to the 2D case in this 
proof (the extension to d > 2 is rather straightforward). There is an evident bound 
11V||u ||^2 < c||zt||x ,2 which implies ||u ||£2 < C\\u \\ L 2 by a Poincare type inequality [ ]. 
To continue, let us change the coordinates on Ll and suppose that there exist new 
coordinates (^ 1 ,^ 2 ) so that D becomes the unit square = (0, l ) 2 and V|| becomes 
cr(^i, ^ 2)^^7 with some positive function a. Problem (46) is written in these new 
coordinates as 

f dv dw f 

/ N WFWF d ^ ld ^ 2= l Juwd £id£ 2 , 

where J = J(£ 1 ,^ 2 ) is the Jacobian and N = N(£ 1 ,^ 2 ) = -A|| J ce 2 , which are positive 
functions given by the geometry. 
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Let us now replace here w by with arbitrary and sufficiently smooth function 
oj such that w = 0at£i = 0 and at 6 = 1. Integration by parts with respect to 6 
yields then 

f d 2 v du 

■ d^dh + / N——— d^di 2 = 

Jn e 0606 06 


f dN dv du: 

ln e 06 06 06 


d(Ju ) 

~d^~ 


ud^\d^2- 


Noting that uj is not differentiated in the last formula wrt £i we can use den¬ 
sity arguments and extend this relation to a broader class of test functions uj, not 
necessarily vanishing at £i = 0,1. In particular, we can now set uj = JjL and get 


in, 


( d 2 v \ 2 

w fefe «= 


f d{Ju) dv 

fe 06 0£i 


dtidt 2 - 


I" ON dv d 2 v 

fe 06 <9£ 2 0 £i0£2 


dtid& . 


This, reminding 




9£ 2 


L 2 


< c||V||n||x ,2 < c| 1 ^ 11 ^ 2 , entails by Young inequality 


d 2 


< cqfellfe + - 

L 2 7 


<9u 


06 


L 2 


T" C| |It| 1^2 


<96 <96 

with a fixed constant c > 0 and arbitrary 7 > 0 . 

Applying a Poincare type inequality to we can write V 6 £ (0,1) 


(47) 


L 


1 f dv 

\06 


d6 < C 


:(&)’*■«= (i?s-,4 ■« 


Remind that w E 7l, which means 

/ fe6,6M6,6fe6 = 0 V6 e (o, 1 ), 

Jo 

or, after differentiation wrt £ 1 , 

fetfefetH- 0 v£i£(o ' i) ' 

Relation (48) can be now rewritten as 


( dv 


- 


V06 


/' 1 / (9 2 t; \ 2 

‘* sc L (kje) ** c ' 


</: 


<9J 

06 


ud6 


which, after integrating over £1 E (0,1), with the aid of (47) and the the bound 
IMIl 2 < C\\u\\ L 2 , gives 




06 


< c 


L 2 


a 2 


0606 


dr 

+ C\\v\\ 2 L2 <C n \\u\\ 2 Hl + — 
L 2 7 


dv 


06 


L 2 




This implies, taking 7 sufficiently big. 

dv 


06 


L 2 


< C||u||#i , 


which gives the desired result since, as already noted, 


dv 




L 2 


< c||rt||x, 2 . 
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Corollary 19 Let £ e A. Then one has |£|y' < c|£|* with some constant c > 0. 


Proof. One can immediately see that |£|y' = |u|y where u E V solves 


(Vu, Vw) = (£,w), Vw€V. 


( 49 ) 


This means in particular that £ = —A u. Let now v € A be the solution 
corresponding to u solution to (49). Lemma 18 implies thus that v € iL 1 (fl) 
has 


l£lv' = Mv = 


(—Art, u) 


\u\ v 


Hm 


a\\{v,0 ^ a \\(^ v ) 

-< c -fp-fr- 

mm IMI m 


< c|«eU- 


to (46), 
and one 


3.2 Numerical analysis for the stabilized AP-scheme 

Let us introduce a mesh Th on 0 consisting of triangles (resp. rectangles) of maximal 
size h and let 14 C V be the space of (resp. Qk) finite elements on 74 We want 
now to discretize the stabilized problem (35) and remark that we can use 14 for both 
variables u and 4 We are thus looking for a discrete solution {u £ ^ a , € 14 x 14 of 


(AP S ) 


£,<T 


a(%’ CT , v h ) + (1 - e)a|| (^ a , v h ) = (/, v h ), Vv h € V h 
a l|( M h <T ’ w h) ~ £a\\{Z,h a ,Wh) ~ o-(^’ a ,w h ) = 0, Vw h € 14 


(50) 


Let us decompose now 14 = G4 © Ah with Gh = 14 n Q = 14 © Qc and Ah 
being the L 2 —orthogonal complement of Gh- Taking test functions from Gh in the 
second equation of (50), we see that £, £ h ’ a € Ah so that this problem can be in fact 
equivalently rewritten as: Find (u £ ^ a ,^ £ h ’ a ) € 14 x Ah such that 


{AP's) 


,/ \£,cr 


a ( u h a i v h) + (1 - 4> a || (^ ,<T , Vh) = (/, v h ) , Vv h € 14 
a\\ w h ) - ea\\ {^ a , w h ) ~ ,w h ) = 0, \/w h € A h - 


(51) 


The advantage of the last reformulation is purely analytical, as we can now reintro¬ 
duce the mesh dependent norm on Ah 


|6iU = sup 
v h ev h 


a\\(€h,Vh) 
\vh\v 


We now equip the space Xh ■= 14 x Ah with the norm \\uh, £h\\x e a h ■■= (\ u h\v + 
\^\lh+^ h \l+^ h \l 2 ) 1 / 2 . By Lemma 12, the bilinear form G C) „ is continuous on 
(X h , ||-, -| \x S)a<h ) and enjoys the inf-sup property 


inf 


(u hl \ h )ex h (, Vh ,wJ)ex h \uh,ih\x e ^ A \vh-,Wh\x l 


C'£,o-{{'Vh, 4i)j {Vhi VJh)) 


> A 


£,(T,h 


(52) 


with a constant /3 > 0 that does not depend neither on the mesh nor on e and u. 
This implies 
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Theorem 20 (Discrete Existence/Uniqueness/a, e-Convergences) The dis¬ 
crete AP-problem (50) admits a unique solution {u/ a ,^ a ) sVjjX Vh, satisfying 


u 


£,a j-£,( 
h >s/i 


< - 


p 


Moreover, for any e > 0 fixed, one has the convergences 


£,(T 


v £>0 V t 

fo—»o u h ; fo—?, 

where («f’°,£f’°) € Vh x Ah is the unique solution to (51) with a 


/-£,0 


0. We also have 


£,0 


“4-X) u 


0,0 


d £,0 

>•h 4—>-0 si 


0,0 

h 


where (u^’ 0 , £jj’°) £ Vh x Ah is the unique solution to (51) with e = a = 0. 

The condition number of the matrix corresponding to problem (50) is bounded by 
a constant that depends on a but not on £ (assuming that the same bases of Vh. and 
Lh are chosen for all values of e,a). 


Proof. The existence and uniqueness of a solution (iq) ,<T , £^ ,<T ) G Vh x Vh for each 
e > 0, cj > 0 is a simple consequence of the BNB theorem [ ]. As mentioned already 
this solution lies in fact in V/, x Ah and it is thus also the solution to (51). By the 
same arguments, the latter problem admits a unique solution {v// 7 ,t/) a ) G Vh x. Ah 
also in the case er = 0. To prove the convergence (v/ 77 , f/ 77 ) —X (u^’ 0 ,^’°) for e > 0 
fixed, we observe 


C E 


// £,CT £,0 /-< 

r(K ~ U h >4 ~i 


‘‘h 

f -£,0 


e,0i 

h > 


, (v h , W h )) = Wh) < crUh 


£,0 | 


< \\LA\ v h,Wh\\x t 


e,cr,h "> 


L 2 \Wh\\L 2 

V(v h ,w h ) G 14 x A h , 


and we conclude using the discrete inf-sup property for C E _ a . The proof of the other 
convergence £->0 while a = 0 is done exactly in the same way as in Theorem 6. 

We turn now to the study of the condition number. We recall the notations from 
the proof of Theorem 6 with the only change that there is no longer the space Lh ., 
which has been replaced by 14. In particular, the constants fi q , v q , //*, u * are now 
evaluated on 14 instead of Lh and one can have fi q = /r* = 0. Denoting by ft and h 
the minimal and maximal eigenvalues of the mass matrix (</>“, (j>j)L 2 (ci) we conclude 
for any e, a > 0 and any = ( u h ,£h) £V h xV h 

m.h/p? u ,a(i 2 )\\m< \\* h \\ 2 x ^ h <max{i/2,i/2 + ei/J + <ri>2}||$|| 2 . 

Introducing the matrix of problem (50), denoted by A £,rr , and repeating the calcula¬ 
tions of Theorem 6 we obtain 


||A £ ’l 2 < M max{^, v 2 + ev 2 + ah 2 } , 
for any $ G M. N , so that one has finally 

cond 2 {A e ' a ) = 11A e,cr 11 2 11(A e, °")— 1 11 2 < 
which is an e-independent bound. ■ 


ft minimafi 2 )\\A>\\ 2 < ||A e ’ CT $|| 2 , 

M max{^, v 2 + ev 2 + ah 2 } 
/3mm(p 2 ,a(i 2 ) 
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Remark 21 As already observed in Remark 7 we have 

fi ~ v ~ h, p u ~ h and n u ~ u q ~ z/* ~ 1. 

Hence, assuming e, cr G [0,1], one obtains 

cond 2 {A £ ' a ) < -j-2 . 

Theorem 22 (h- Convergence) Let k > 1 and If &e t/ie Pk or Qk finite element 
space on a regular mesh Th■ Suppose moreover that problem (35) has the solution 
u £,a £ H k+1 (Ll), f ,or € lf +1 (f2) and problem ( 40 ) has a solution u £ € idffl), 
f € HThen 

\u £ - u e f \ H i < Ch k (\u e,c \ H k+i + |f ■‘V+i) + Ca|f | H i , (53) 

with a constant C > 0 that depends neither on the mesh, nor on £ or a. 

Proof. In the same way as in the inflow case we prove that 

\u"' a — uf |j^i < Ch k {\u £ + If |j? fe +i)- 

It remains to invoke Lemma 14 and the triangle inequality to conclude. ■ 

Remark 23 The error estimate (53) would be of course useless if the norms \u e,a \ H k+i, 
\ffi' a \ H k+i were dependent on e and a. Fortunately, it is not the case. We expect in¬ 
deed that |n e,CT |^fc+i is bounded uniformly in e by the norm of f in H k ^ 1 (Ll) and 
|f ,a \ H k+i is bounded uniformly in e by the norm of f in H k+1 (Ll). This can be easily 
proved in the case of simple aligned geometry, see Appendix A. 

Remark 24 One can also easily obtain 

+ <r\ ri H .] 

which degenerates as in the inflow case, as e goes to 0. Again, we are not sure, if 
this estimate is sharp, but we recall that f ,<T is an auxiliary variable, without any 
intrinsic meaning. 


4 Numerical tests 


Let us now study numerically both AP-reformulations, the inflow as well as the 
stabilized one. We consider in the following a square computational domain Q = 
[0,1] x [0,1] and the non-uniform and not coordinate-aligned b field: 


B=( a(2y — 1) cos(-7rx) + ir \ 

B | ’ \ na(y 2 — y) sin(7rx) ) ’ 


(54) 
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as well as a sample function u° which is constant in the direction of the 6-field: 

u° = sin (iry + a(y 2 — y) cos(7rx)) . (55) 

Here a > 0 is a parameter to be fixed in the following different test cases and describes 
the variations of 6. We choose u° to be the e —>• 0 limit solution of the anisotropic 
problem ( P) £ , hence solution of (12), and construct an exact solution of (10) by 
adding a perturbation proportional to e, i.e. 

u £ = sin (iry + a(y 2 — y) cos(7rx)) + e cos (27 tx) sin (ny + a(y 2 — y) cos(7rx)) . (56) 
Note that the auxiliary variable q £ , solution of (15), is in this case equal to 

q £ = cos ( 27 tx) sin (7 xy + a(y 2 — y) cos(7rx)) — sin (7 xy + a(y 2 — y) cos(7rx)) . (57) 

Finally, we compute the right hand side accordingly and have thus constructed 
an exact solution of problem (10). All simulations (unless stated otherwise) were 
performed using a Q 2 finite element method. 

Aim of this section is to study and validate from a numerical point of view the error 
estimates established in the last two sections. In particular, we investigate firstly the 
error introduced by the stabilization procedure in the ( APs) £ ’ a formulation, meaning 
the cr-convergence estimate of (45) in Theorem 14 is verified numerically. Then the h- 
convergence of both methods is studied and the estimates (32) and (53) are confirmed 
in both anisotropic (e <C 1) and isotropic (e ~ 1) regimes. Next, we show that both 
methods are Asymptotic-Preserving in the parameter e. The conditioning of the 
corresponding linear systems appear effectively to scale in agreement with Remarks 
7 and 21. Finally, the case of a less regular force term /, belonging merely to L 2 (Q) 
(and not to i7 1 (H)) is studied — the convergence of the schemes is tested beyond 
the validity of Theorems 8 and 22. 

4.1 Stabilization error (e > 0, h > 0 fixed, a -A- 0) 

Let us start by studying the error introduced by a stabilization term proportional 
to a in the ( APs) £ ' a reformulation, in particular we shall estimate numerically for 
fixed s > 0 and h > 0 the L 2 - resp. H l -e rrors between the exact solution u £ con¬ 
structed in (56) and the numerical stabilized solution u e ^ a , solution of (35) or (50), 
i.e. | \ u £ — u £ h a \\. The mesh size h is fixed to 0.01. Numerical simulations are per¬ 
formed for the stabilization constant a varying from 1 to 10~ 15 , considering three 
different regimes : no anisotropy (s = 1), strong anisotropy with direction aligned 
with the coordinate system (e = 10 -10 , a = 0) and strong anisotropy with variable 
direction (e = 10~ 10 , a = 2). The L 2 - and 77 1 -errors are presented as a function of 
a in Figure 3. 

In the first regime, with no anisotropy present in the system, the stabilization con¬ 
stant does not influence the precision at all. Indeed, it is exactly what is expected as 
the terms involving £ £,a do not appear for e = 1 in the first equation of ( APs) £ ' a . In 
the second regime, with strong and aligned anisotropy, the precision of the scheme is 
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Figure 3: Absolute error \\u £ — u^Wl 2 (on the left) and \\u £ — u^Wh' (on the right) 
with respect to the exact solution u £ , as a function of cr, for h = 0.01 and three regimes 
: no anisotropy, strong and aligned anisotropy as well as strong anisotropy with variable 
direction. 

influenced by the stabilization procedure only for a -values greater than 10 -5 in the 
L 2 -norrn and greater than 10 -3 in the H 1 -norm. The error dependence in a is here 
linear, according to the Theorem 14, and is explained simply by the fact that the 
stabilization influences the results for large a > 0. 

For e-values smaller than these critical values the accuracy of the scheme in both 
norms remains unchanged and is given only by the mesh size. This holds true even if 
the value of the stabilization constant is close to the machine precision (10~ 15 ) and 
can be explained by the fact that we are in an aligned test-case. Indeed, normally 
for <7 —>• 0 the error should increase, due to the non-uniqueness of the ^-solution. 

Here we are however plotting the error corresponding to the rh-function, which is 
uniquely determined. The non-uniqueness of £ £,<T steps in only in the not-aligned 
case, which is our third regime of strong anisotropy with variable direction. In this 
case, the curves show an expected cr-behaviour, the optimal value of cr being between 
10 -8 and 1CH 5 for the L 2 -error and between 10 _1 ° and 10~ 3 for the //'-norm. To 
explain this, observe that in the limit a —» 0 the auxiliary function £ £ ' a is uniquely 
determined up to a constant on the field lines. This constant will normally not in¬ 
terfere in the computation of u £,a , as only the parallel derivatives of £ £,a are present 
in the u £,CT -equation. However, if the mesh is not aligned with the field lines, this 
parallel derivative mixes the directions, introducing errors which lead to the observed 
behaviour of the error as cr —r 0. 

Having tested the cr-dependence of the error ||u £ — u^ a || for fixed h > 0, we are 
now interested in how these curves remodel for different /t-meshes. The u-convergence 
is hence compared for different mesh sizes in the most difficult setting, that is to say 
when a strong anisotropy with variable direction is present in the system (e = 10~ 10 , 
a = 2). Numerical simulations were performed for the mesh size ranging from 0.1 to 
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(a) L 2 error 


(b) H 1 error 


Figure 4: Absolute error \\u £ — (on the left) and \\u £ — u^Wh 1 (on the right) as a 

function of a, for different values of h and for £ = 10 -10 and a = 2. 


0.003125. Cumulative results are presented on the Figure 4. The plateau for which 
the accuracy of the scheme does not depend on the stabilization parameter is clearly 
dependent on the mesh size. As a consequence, the value of a should be clearly made 
mesh dependent. We observe that in the case of Q 2 finite elements the upper and 
lower bounds for the optimal value scale like h 3 and h 4 for the L 2 -error, while for the 
H l -e rror the respective scaling is approximately h 2 and /i 6 . It is therefore reasonable 
to put a = h 3 (or a = h 2 if one is interested in the H l - precision only). Note that 
this scaling depends on the finite element method used. In general, if a (or Qk) 
method is used, the optimal choice of a is h k+1 , which ensures optimal h-convergence 
of the method in the L 2 -norm. 

4.2 /^-convergence (e > 0 fixed, a = h 3 , h —> 0) 

Let us now turn our attention to the /i-convergence of both Asymptotic Preserving 
reformulations (21) and (35). Since the AP-scheme with inflow boundary conditions 
was studied in a previous work [ ] we are mainly interested in the behaviour of 
the scheme with stabilization. As in the previous subsection, numerical tests are 
performed in three regimes : an isotropic one (e = 1 and a = 0) and two anisotropic 
regimes (e = 10~ 10 and a = 0 or a = 2 ). The stabilization coefficient is set to 
a = h 3 as a consequence of the last subsection. The convergence rate in the L 2 - and 
fL 1 -norms is presented on Figure 5. As expected the optimal convergence rate (of a 
Q 2 -FEM) is found in both norms. Next we compare the results with the convergence 
of the AP-scheme with inflow boundary conditions in Tables 1 and 2. Note that in 
the case of no anisotropy or anisotropy aligned with the coordinate system (a = 0), 
both schemes give quasi exactly the same precision for both L 2 and fL 1 -norms. In 
the last regime the stabilized scheme is slightly less accurate compared to the ( APi n ) £ 
scheme. A small loss of the convergence rate of the stabilized scheme is observed for 
the smallest mesh size in both norms. 
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h 

e = 1, a = 0 

£ = 10~ 10 , a = 0 

£ = IQ' 10 , a = 2 

(■ AP in ) £ 

{AP s ) £ ’° 

(■ AP m ) £ 

(. AP s ) £ ' a 

(AP in ) £ 

(AP S ) £ ' C 

0.1 

5.39 x 10" 3 

5.39 x 10" 3 

1.19 x 10" 3 

1.19 x 10" 3 

2.81 x 10" 3 

2.18 x 10" 3 

0.05 

6.97 x 10~ 4 

6.97 x 10 -4 

1.49 x 10” 4 

1.49 x 10 -4 

3.16 x 10 -4 

2.87 x 10 -4 

0.025 

8.79 x 10~ 5 

8.79 x 10" 5 

1.86 x 10" 5 

1.86 x 10" 5 

3.77 x 10" 5 

3.53 x 10" 5 

0.0125 

1.10 x 10" 5 

1.10 x 10" 5 

2.33 x 10~ 6 

2.33 x 10" 6 

4.57 x 10" 6 

4.31 x 10" 6 

0.00625 

1.38 x 10” 6 

1.38 x 10" 6 

2.91 x 10" 7 

2.91 x 10" 7 

5.60 x 10" 7 

5.29 x 10" 7 

0.003125 

1.72 x HT 7 

1.72 x 10" 7 

3.64 x 10" 8 

3.64 x 10" 8 

6.89 x 10" 8 

6.52 x 10" 8 

0.0015625 

2.15 x 10" 8 

2.15 x 10" 8 

5.51 x 10~ 9 

4.78 x 10" 9 

1.07 x 10" 9 

8.05 x 10" 9 


Table 1: Comparison of the L 2 relative precision \\u £ — u] II^/II^^Hl 2 of both reformu¬ 
lations in both isotropic and anisotropic regimes for different mesh sizes and stabilization 
constant set to a = h 3 . 


h 

e = 1, a = 0 

e = 10 ” 10 , a = 0 

e = 10 " 10 , a = 2 

(. AP in y 

(■ AP S ) £ ' C 

(. AP in y 

(■ AP s y a 

(. AP in y 

(■ AP S ) £ ’ ff 

0.1 

4.48 x 10" 2 

4.48 x 10" 2 

1.46 x 10" 2 

1.46 x 10" 2 

2.44 x 10" 2 

2.33 x 10" 2 

0.05 

1.13 x 10” 2 

1.13 x 10" 2 

3.67 x 10" 3 

3.67 x 10" 3 

6.34 x 10" 3 

6.12 x 10" 3 

0.025 

2.84 x 10" 3 

2.84 x 10" 3 

9.19 x 10~ 4 

9.19 x 10" 4 

1.60 x HT 3 

1.54 x 10" 3 

0.0125 

7.11 x 10~ 4 

7.11 x 10" 4 

2.30 x 10~ 4 

2.30 x 10" 4 

3.99 x 10" 4 

3.83 x 10" 4 

0.00625 

1.78 x 10 -4 

1.78 x 10 -4 

5.75 x 10" 5 

5.75 x 10" 5 

9.93 x 10" 5 

9.53 x 10" 5 

0.003125 

4.45 x 10" 5 

4.45 x 10" 5 

1.44 x 10" 5 

1.44 x 10" 5 

2.46 x 10" 5 

2.37 x 10" 5 

0.0015625 

1.11 x 10" 5 

1.11 x 10" 5 

3.76 x 10" 6 

3.76 x 10" 6 

6.08 x 10" 6 

5.87 x 10" 6 


Table 2: Comparison of the H 1 relative precision || u £ — ^’ <T ||ff 1 /||^’ <T ||ir 1 of both reformu¬ 
lations in both isotropic and anisotropic regimes for different mesh sizes and stabilization 
constant set to a = h 3 . 
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(a) L 2 error 



Figure 5: Absolute error \\u £ — % ,<t ||l 2 (on the left) and \\u £ — (on the right) as a 

function of h and fixed a = h 3 . One isotropic regime (e = 1, a — 0) and two anisotropic 
ones: e = 10“ 10 and a = 0 or a = 2 are investigated. The optimal convergence rate is 
found. 


4.3 AP-property (h > 0 fixed, a = h 3 , £ -A 0) 

Next, we test if both schemes are indeed Asymptotic Preserving as e —>• 0. The mesh 
size is fixed to h = 0.01, a is set to a = /i 3 and numerical simulations are performed 
for a variable anisotropy direction (a = 2) with an anisotropy strength e varying 
from 10“ 20 to 10. Both schemes exhibit the desired property, as shown in Figure 
6. In particular, the absolute error for both reformulations and in both norms is 
independent of e (for £ < 0.1). The error curves are practically indistinguishable. 
For large e- values, the errors are increasing due to the fact that the here presented 
schemes are designed to cope with e <C 1 singularities. 


4.4 Matrix conditioning 

Finally, let us now turn our attention to the conditioning of the matrices associated 
with the numerical resolution of both schemes (APi n ) £ resp. (APs) £,a . The strong 
anisotropy case with variable direction (a = 2) is considered for different mesh sizes 
h > 0. The stabilization constant a is set to h 3 in the ( APg) £ ' a reformulation and 
the anisotropy strength e is set to 10“ 10 . Sparse matrices were assembled in every 
case and the condition number was estimated using the matlab function condestQ 
returning the estimate of corul\. The results are displayed on Figure 7. As expected, 
the conditioning scales as 1/h 4 for the inflow reformulation and as 1/ah 2 = 1/h 5 
for the stabilized method. The first method results in better conditioned matrices in 
this setting. However, if one is interested mainly in the H 1 precision the stabilization 
constant a could be set to h? resulting in a conditioning proportional to 1 / h 4 for the 
stabilized method, discretized with the (Qb finite elements. 
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(b) H 1 error 


Figure 6: Absolute error | u £ ex — u^ m \ \ jy. (on the left) and | \u e ex — u^ m \ \#i (on the right) as 
a function of e for an anisotropy not aligned with the coordinate system (a = 0) and the 
mesh size h = 0.01, a = h 3 . The error curves are superposed, both schemes show similar 
accuracy independently of e. 



Figure 7: Conditioning ( cond \) of the matrices associated with both AP schemes as a 
function of the mesh size for strong and nonaligned anisotropy (e = 10 -10 , a = 2). The 
predicted scaling is found. 
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4.5 The case of f ^ H 1 (S1) 

Aim of this subsection is to investigate the error estimates in a case where the right 
hand side / is less regular than supposed in the theoretical part of the last two 
sections. All simulations in this section are preformed with a Qi finite element 
method and the stabilization parameter in the (APs) £,<7 -formulation is set to h 2 . In 
this case we have the h-estimates (32) resp. (53) with k = 1 and we recall Remark 9 
resp. 23. Let us now choose u° to be defined by 

u° = ( {y + a(y 2 — y) cos(7rx)/7r) 2 In (y + a{y 2 — y) cos(7rx)/7r) — 1.5 J 

^ * (58) 

+ 7.5 (y + a(y 2 - y) cos(7rx)/7r) , 

so that the right hand side for the limit problem is a function that belongs to L 2 
and not to H 1 . If a = 0 (the field is aligned), then the right hand side of the limit 
problem equals to In y. 

We remind that in view of our theoretical result, the iT-norms of q £ and £ £,CT are 
not guaranteed to be bounded if the force term is not H l {0). Nothing can be said 
on the convergence of the numerical methods in this test case since the right hand 
side (58) is not in H 2 ( fl). We consider two anisotropic regimes (e = 10 -10 ): with 
anisotropy direction aligned with the coordinate system (a = 0) resp. with variable 
direction (a = 2). 

Numerical simulations show that the H 1 -norms of q | and £^ ,<T grow as h ap¬ 
proaches 0 for the variable anisotropy direction. This seems to confirm our expec¬ 
tations. On the other hand, the H 1 -norm remains constant when the anisotropy is 
aligned with the coordinate system. The L 2 -norm seems to be bounded regardless of 
the method for both studied regimes. The results are displayed on Figure 8. To our 
surprise, optimal convergence rate of u e h and u^ £ is conserved in the tested h range 
— see Figure 9. 


5 Conclusions 

A detailed numerical analysis of some asymptotic-preserving numerical schemes, de¬ 
signed to cope with highly anisotropic elliptic problems, was carried out in the present 
work. In particular, we have shown rigorously that in the limit regimes where tra¬ 
ditional schemes become inadequate, AP-schemes are perfectly able to capture the 
macroscopic behavior of the solution. Convergence results for the schemes were 
proven, with an accuracy and stability which are shown to be e-independent, e being 
the perturbation parameter responsible for the stiffness of the problem. The devel¬ 
opment of AP-schemes is based on asymptotic arguments and permit hence to create 
a link between the various scales in the considered problem, while the numerical 
parameters remain independent on the stiffness parameter. 
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(b) H 1 norm 


Figure 8: L 2 (on the left) and H 1 (on the right) norms of £ e,a and q £ as a function of the 
mesh size h > 0, for £ = ICO 10 , a = h 2 and a = 0 or a = 2. The H 1 -norms of both 
auxiliary variables increase with decreasing mesh size for variable direction of anisotropy. 
The L 2 -norms seem to be bounded. 




(a) L 2 norm (b) H 1 norm 


Figure 9: Absolute error \\u £ — u e ^ a \\ L 2 (on the left) and \\u £ — (on the right) as 

a function of the mesh size h > 0, for £ = 10 10 , a = h 2 and a = 0 or a = 2. Optimal 
convergence rate is observed for both methods and both anisotropy configurations. 
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A The regularity of the solution in the case of 
a simple geometry 

Consider the case of 12 = (0,7r) 2 and the field b looking upwards, i.e. b = e y . More¬ 
over let .A|| = 1 and A± = Id. We want to explore in this Appendix the regularity of 
the solution (u £,a ,£ £,CT ) to (35) when / belongs to H S (Q) and considering an aligned 
geometry case. 

To start, let us first remark that the functions {^/2pn sin kx}k> 1 as well as {y / 2/7r cos lx}i> o 
form an orthogonal basis in L 2 (0, vr) [2], such that each / G L 2 (Q) can now be written 
under the form 


OO OO 

f(x,y) = EE fkl sin kx cos ly , { fkl}k,le N C l 2 , 

fe=l 1=0 


which implies immediately that 

CO CO J 

u £ ’ a (x, y) = ^2^2 — --- OF f kl sin kx cos l V ’ 

k =i i=o k2 + l2 + 1FT7 

oo co j2 

t) = E E ( £i 2 + g )(t 2 +i 2 ) + (i_ e)i 4 /» sin kx cos ly. 
Now, if / € £P(f2), Parseval’s equality permits to show that 

CO OO 

i/i^~ES > 2 + * 2 r/* 2 z, 

k= iz=o 


so that 


| a |//s+2 


(A 2 + / 


2\s+2 






fc=lz=o 


I/I 


2 

// 3 > 


If’' 7 1 2 


H 3 


EE 


(A: 2 + l 2 ) s l 4 


k = i /=i (( e ^ 2 + cr )(^ 2 + l 2 ) + (1 — e )^ 4 ) 
Moreover, in the case cr = 0 one has 


2 fki < EE( fc2+ i 2 yfki ~ i/i^ s • 


fc=lZ=1 


|(9 ra u 


ei2 

I// 3 


(V.2 I ]2\si4 00 00 

EE ; A fi, <EE(‘ ! +~^ 2 i/i«- 


fc=l z=i ^fc 2 + ^ 


fc=l 1=0 


In conclusion, if f G H S (Q) then u £,rT G H S+2 (Q), f ,f7 € H S (Q) and d yy u £ G H S (Q) 
and there is a constant C > 0 independent of e and a such that 


I U" ,CT \h s +2 


<C\f\ 


H 3 


|£ £ ’°1h 3 < C\f\ H * and \d yy u e \ H 3 < Ce\f\ H 3 ■ 
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The same estimates hold true in the inflow case, problem (15), i.e. when q £ is 
associated with zero boundary conditions on the inflow part. Indeed, the link between 
q £ and £ £, ° can be explicited as 


q £ (x,y)=^°(x,y)-^°(x, 0). 


Hence, it suffices to study the regularity of the trace function x £ ( x ) := £ £ ’°(x, 0). We 
have 

\ \ ^2 

= E E e p (t 2 + |2 ) + (1 _ £ ),4 /« Si “ kx . 

implying 


lx £|2 


H°(T in ) 


E fc2s E 


fkl 


k =1 

oo 


“ el 2 (k 2 + l 2 ) + (1 — e)l A 


/«) ^E fc2 1E |2 


fc=i \i=i 


< E^(Ef](E/^I^EE(^«h/i^ 


fc=i 


/=! 


/=! 


k =1 2=1 


We conclude thus x £ € -H s (rj n ) so that g e € H s {£i) with the e —independent estimate 
W £ \h s < C\f\n s - 


B On the discrete inf-sup condition 


As mentioned earlier, the numerical analysis in this paper would be more convenient, 
if the discrete inf-sup condition (29) were true, i.e. 


. , aJq h ,v h ) 

mt sup -j—j—j—j— 

qh&L hVh ev h \qh\*\Vh\v 


> a , 


with a mesh independent a > 0. In more explicit form, this means 


w r a»(q h ,v h ) a\\(q h ,v) 

yqh € L h : sup —:—:- > a sup ——:- 

v h ev h \Vh\v vev \v\v 


(59) 


(60) 


We show first that (59) holds true in a simple aligned geometry. Secondly, we 
provide a numerical study of a non-aligned case where (59) turns out to be false. 


B.l The case of the aligned geometry 

Assume H = (0, L x ) x (0 ,L y ) and b = e<i- Choose 14 as the finite element space Qk 
(k > 1) on a rectangular grid Th aligned with the coordinate axes. More precisely, 
we choose some node points 0 = .To < x\ < ■ ■ ■ < xn x = L x on [0, L x \ and 0 = 
Ho < Vi < • • • < ljN y = L y on [0, L y ] with all the steps of order h and introduce Th 
as the collection of rectangles [xi,Xi+ 1 ] x [yj,yj+ 1 ] that constitutes a partition of O. 
Moreover, £h shall denote the set of all the edges of the mesh Th with the exception 
of those lying on To. 
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Our strategy to prove (59) is to use Verfiirth’s trick [14] by first establishing the 
inf-sup conditions with respect to an auxiliary mesh dependent norm and then going 
to the original norm with the aid of Clement interpolation. Thus, we want to prove 
first that for all qh E Lh there exists Vh £ Vh such that 

Nh(qh) := (l £ X IM|2<fa+ L |iW|2<fa ) £ c ~\ wt <61) 

where [d y qh] denotes the jump of d y qh across an edge E E if the edge is internal, 
and [d y qh\=d y qh on an edge E lying on the boundary T. 

A convenient reformulation of this is: For all qy t E L h there exist Vh E Vh such 
that 

\\vh \\ 2 L 2 < C x Nl{q h ) and a\\{q h ,v h ) > C 2 N%(q h ). (62) 

The quantities above can be written in a more explicit manner as 


N 2 h (q h ) 

I kills 

a\\{qh,Vh) 


1 ' ' v . ‘' y rx x rm 

I [d y q h }\ 2 (yi)dx + Y1 / / I dy y q h \ 2 dydx, 

h i=0 Jo i=0 Jo J y*-1 

r L x r Ly 

/ / v\dydx , 

Jo Jo 


y rL x ryi 


N , 


v r L x 


N , 


kkk yi)vh{x-, Vi)dx - ^ 


^ pL x ryi 


d yy q h v h dydx. 


i =0 


i=i J0 J yi-i 


The construction of Vh is particularly easy in the case of bilinear finite elements 
(k = 1): we can take Vh E 14 such that for all i = 0,..., N y 


Vh(x, yi) = -[d y q h ](x,yi) . 


(63) 


Then the second inequality in (62) becomes equality with C 2 = 1. Moreover, one 
easily gets by a scaling argument 


rL y N v 

/ Vh( x ,y)dy < C 1/1 Vk)k,2/i) 

i=Q 


(64) 


for all x E [0, L x \ which gives the first inequality in (62). 

We describe now a more complicated construction in the case k > 2. For a given 
qh € Lh, we construct Vh E Vh as v h = v y ' > + where E Vh is defined for any 
x E [0, L x \ and any y E [y,;_ 1 , i = 1,..., JVj, by 

= -d yy qh(x,y )———— 


36 




( 2 ) 

and v h &Vh is such that 

v h\ x ,yi) = j L [dy<lh\(x,yi), i = 0,...,N x , xe[0,L x ] 

(2) 

and v y h \yG[yi-i,yi] f° r any x fixed is a polynomial of degree k that is orthogonal in 
L 2 (yi_i,yi) to all the polynomials of degree ^ k — 2. We get for this Vh 


N, 


v rL x 


N, 


(Qh,v h ) = 7^ / \[ d yQh]\ 2 {yi)dx+^2 


y pE x ryi 


h ‘ ” / n 

i =0 


i= o Jo J yt-i 


\d yy qh\ 


2(y-yi-i)(yi ~v) 

h 2 


dydx 


This yields immediately the second inequality in (62) with some C 2 > 0. In order to 
prove the first inequality in (62) we employ again a scaling inequality of type (64): 
for any x E [0, L x ] 

rLy pEy pEy I pVi 1 

/ v\dy< 2/ \v^\ 2 dy+2 \v^ ] \ 2 dy < C x y / \d yy q h \ 2 dy + - y[^,©J 2 (j/i) 

Jo Jo Jo \ i=1 J yi -1 n i=0 

Now, (61) being established, take any qh € Lh■ By the definition of the norm 
| • I*, there exists v E V such that |u|y = \qh\* and a\\(qh,v) = \qh\ 2 - Let Vh E 14 be 
the Clement interpolant of v such that [7] 


\v — 


*yz,2(fi) < Ch\v Iv, \vh\v < C'l^lv, ( y ll u ^IIl 2 (E) ) - ^|u| v . 

E^Eh 


Observe that 

\Qh\l = a\\(q h ,v) = a\\(q h ,v-v h ) + a\\(q h ,v h ) 

f ~ ^ 1 a \\{QhiVh) ( ~ \ X r>\ \ a \\{Qh,Vh) 

< a\\{q h ,v-v h ) + \v h \v sup — :—:- <a\\(q h ,v-v h ) + G\q h \* sup — : — :- 

v h ev h \vh\v v h ev h K|y 

Integrating by parts element by element in the first term of the last line yields 


a\\(q h ,v-v h ) = V' / Vug/, • V\\(v - v h )dx 
777 Jk 


KeT h ' 


y / h|.-V|| q h ](y-v h )ds- V / (V 11 - V M q h ){y-v h )dx 
77c Je 777 Jk 


E^Sh ' 


KeT h ' 


z aw \^r>u\\ a \\(9h,v h ) a\\(q h ,v h ) 

< Nh(qh) < G/r|r;|y sup ——-- < G|u|y sup 


IKIIl 2 


^ "JC I I 

v h ev h \Vh Iv 


We have used here (61) and the standard inverse inequality. This enables us to 
conclude 

, ,2/w 1 “II (<lh,v h ) , a\\(q h ,v h ) a\\{q h ,v h ) 

\qh\* < C\V\v sup -:-;-b C\qh\* sup -:-:- < C\q h \* sup -;-;- , 

v h ev h \Vh\v v h (EV h \Vh\v V h &v h \Vh\v 


since qh E 14 and |u|y = \qh\*- The last inequality gives the desired result (59). 
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B.2 A numerical study in the case of a general geometry 

The aim of this section is to investigate the validity of (59) or equivalently (60) in 
a more general context by a series of numerical experiments. As in Appendix B.l, 
we shall assume that = (0 ,L X ) x (0 ,L y ), b = e 2 , however this time the grid is no 
more aligned with the field lines of b. Indeed, we are using here a regular grid made 
of triangles such that their hypotenuses are no longer aligned with b. Numerical 
simulations are performed with FreeFem++ [10]. 

Let Vh be the V\ finite element space on a mesh described above of size h > 0. 
Observe that the first supremum in (60) is attained on v^ € 14 that satisfies 


a(v* h ,w h ) = a\\(q h ,w h ), Vw h <EV h . (65) 

In order to explore the second supremum in (60), we use the finer finite element 

space V ^/21 constructed via V 2 finite elements on mesh of size h/2, i.e. a two-times 

refinement of the mesh above. The goal in introducing this finer space V^ 2 is to 

approximate the infinite-dimensional space in (60). 

* f f 

Consider v^ 2 € V^ 2 that satisfies 

a «/2’ w h/2 ) = a ll (9ft/ 2, w f h/2 ) , Vu { /2 € Vl j2 . (66) 


If (60) holds true, than we have 


V5/1 € Lh : 


J h lv 


* f 


V h/2 

V 


> a. 


Unfortunately, this is false as shown in the following numerical experiment. Let 
L x = L y = 1 and let us choose on each mesh of size h = l/n the function € 14 
defined by its values at the mess nodes as 


qh(xi,Vj) = XiSin(irnyj/2) , (67) 

where Xi = ih, rjj = jh, i,j = 0, Note that this function satisfies all the 

boundary conditions provided n is even. In Fig. 10 we plot the quantity 

K/2 lv 

computed for such a qh on a series of meshes versus /t = A It shows clearly that the 
constant a in (59) is mesh dependent, i.e. it tends to 0 (in general) when the mesh 
size tends to 0. 
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Figure 10: The ratio } v f\f 

IVI v 


computed for q h given by (67). 
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